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ABSTRACT 



i , 



A dynamic model of a multicomponent-multistage distil- 
lation column with variable holdup and boundary lags is developed, 
and the programming for digital simulation is presented. The open 
loop transient responses for the mcdel undergoing feed composition 
disturbances and reflux change* are studied. These responses are 
then used for synthesizing transfer functions and closing the loop 
for digital control of the overhead product using a variable reflux 
rate. 
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= Dimensionless time = (F/W ) 9 

Lr 

= Dimensionless time increment for integration 




dh . 

J 

(W Lr /F)d T 



= Time rate of change of enthalpy of liquid 
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= Time rate of change of mole fraction of 
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1. Introduction. 



Consideration is given here to the dynamic performance of distil- 
lation columns with the ultimate object being column control studies 
by digital techniques. This is accomplished by mathematically model- 
ing a multicomponent-multistage column and using the model to generate 
column response. The column response is then used as an error source 
so that a control algorithm can generate a corrective action to move 
the system back to its desired control point. The model performance, 
error generation and control action are all done in a digital com- 
puter so that in essence direct digital control is being exercised. 

To have realistic control action, model responses should dupli- 
cate actual plant responses. In order to do this, it is necessary 
to include at the very least column hydraulics and holdup and t delay 
effects in the column end conditions. In this way, variable stage 
holdup is included in the model and the model will respond to changes 
in reflux flow which is one of the possible manipulable variables in 
this work. The work of Peiser and Grover [1] puts forth a model 
similar to the one in this work. 

It is possible to generate model response without these secondary 
effects. Models of this type involve only mass transfer effects but 
can be quite complex if stage mixing and stage efficiency are incorpo- 
rated. In addition it is possible to have interaction between hydro- 
dynamics and mixing and efficiency effects. In the work presented here 
it is assumed that the stages are perfectly mixed and that the ef- 
ficiency is 100%. This in no way limits the model developed since 
allowance for both these effects could easily be incorporated into the 
model . 

Since a model will be of primary importance for response gener- 
ation, it is important to consider all possible physical processes 
which can effect a system response. Unfortunately, a model to do 
this becomes so mathematically unwieldy as to be useless. It has 
become the custom to make simplifying assumptions which produce a 
model that is relatively easy to work with. The sacrifice made is 
loss of realism in responses produced, but this is not too important 
for purposes of study of general system behavior. However, to study 
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system control using models for response generation, it becomes 
necessary to put back in a model some of the "realism" lost by making 
the "usual simplifying assumptions." The problem of what to include 
in a model becomes a difficult one. 

Rosenbrock [2] in his excellent paper on distillation column 
control points out that the difference between observed and generated 
system response at small values of elapsed time is mainly due to 
secondary effects often not included in dynamic models. These second- 
ary effects are the hydrodynamic delays within a column and the delay 
and holdup effects in condenser and reboiler operation and external 
mass transportation. Unfortunately, it is at low values of elapsed 
time that errors for control action may be generated and so a model 
to be used for control studies should include these secondary effects. 

Since distillation has always been a much used operation in 
chemical engineering, the literature is voluminous regarding steady 
state analysis of multicomponent systems and dynamic analysis of 
binary systems. The literature on dynamic analysis of multicomponent 
systems is somewhat less voluminous. In particular, work on multi- 
component systems wherein secondary effects are included is quite 
limited. Some recent publications in this area are the work of Sar- 
gent [3], that of Peiser and Grover [1], that of Moczek et al. [4] and 
finally that of Anisimov. [5] There are many other papers (see for 
instance [6], [7], [8]) besides these, but these are representative. 

As previously mentioned, then, the work presented here is the 
development of a multistage-multicomponent model which allows for 
variable stage holdup, takes account of liquid hydraulics and, 
finally, allows for delays and holdup in the condenser and reboiler. 
This model is used to generate responses to disturbances in the feed 
composition and the reflux flow. The responses are used in two ways 
as follows (1) to abstract the transfer function at any stage or the 
end of the column and (2) to generate the error needed for control 
action by means of a numerical algorithm. The transfer functions 
may be used for analog computer studies cf column control. 

One of the objectives of this work was to see if a numerical 
approach to column control studies could be done in reasonable time 
using digital computers. As might ba expscted, the answer to this 
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is a function of the size of the system studied and the type of 
control involved^ but indications are that work of this type can be 
done in reasonable time at reasonable cost. Another objective was 
to produce as accurate a model as possible with the criterion being 
duplication of actual plant response. Unfortunately, actual plant 
data was not available so that model checking could be done. Never- 
theless) the model as presented and used is accurate in the sense 
that at all times during a transient, the mole fractions sum to unity 
for relatively large values of the independent variable increment. 
This latter point is considered in some detail in the work of Mah 
et al. [9] and Sargent. [3] 
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2. Development of Model. 



General Model 




The above j — stage (j is the stage position index increasing 
vertically) shows all streams crossing the stage boundaries. All 
these streams with the exception of Q carry mass and energy. Q, the 
heat exchanged with the surroundings, is an energy stream only. 

A general dynamic mass balance for component "i" at stage "j" is 

■St 'Wj + Tb <Vi>J = (EinpUt - &UtpUt) mass < L > 

£input = [(Vy,)^ + (LX.). +1 ] + [f^ + 

+ ^ SSFyi VSSF +LsSF \sSF3 

■ K>, + '“i’J * J 

A general dynamic energy balance for the same stage is 

TZ (W h) + (W H) = (Zinput - ^output) (2) 

da v L j da 1 V j v v ' energy v 
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Einput + ^ Lh )j + l] + tv^V + F L H Fl] + tsSF^/SSF + ^Sf'YsSf] 

EoutpuL . [(V H)j + (L h)j ] + Q + [v ssp H vssp + L sgp h LSSp ] 



There will be (M*N + 2) equations of type (1) and (N + 2) equa- 
tions of type (2). The +2 accounts for the boundary condition equa- 
tions. In addition, there will be (M«N + 2) vapor-liquid equilibri- 
um relationships of the form 



ij 



'lAj 



and (N + 2) mole fraction constraint equations 



? y ij = 1 



Pu = 1 



(3) 

(4) 



The equations as stated imply that (a) stage holdups, W and 

W . , are perfectly mixed, (b) equilibrium between vapor and liquid is 
L j 

instantaneously realized at every part of a perfectly mixed stage, 

(c) holdup of vapor and liquid between stages is negligible, i.e. 
downcomer effects are not included, and (d) time to attain fluid 
equilibrium is small compared with that for mass transfer. 

Specific Model 

The general model, while reasonably complete, should be modified 
to reflect only the amount of "realism" required to accomplish the 
purpose of the model. Otherwise, the equations can become unwieldy 
and needlessly complex. With this in mind, the following additional 
assumptions are made about column operation: 

(a) the column operates adiabatically. (Q = 0) 

(b) vapor holdup is negligible compared to liquid holdup up to 
moderate pressures and can be neglected. 

(c) interaction between pressure drop changes, liquid hydraulics 
and mass transfer efficiencies can be nqglected. 

These assumptions can be relaxed and the model modified to in- 
corporate any of these effects should it be desirable to do so. 

Equations (1) and (2) are next rendered dimensionless by dividing 
through by the feed rate, F. A dimensionless time unit, t, is 



15 



introduced, defined by the equation 



T = (F/W ,) 0 



( 5 ) 



where the time, measured in units of tau, is equivalent to the number 

of tray turnovers based on the column feed rate. See Appendix VI for 

t, 0 conversion factor. Further, since it is possible for W . to be 

1*J 

different at every stage and vary with time, a particular W . must be 

h j 

chosen for use with the definition equation (5). A reference stage 

is arbitrarily chosen and the holdup value on that stage at some point 

in time is designated (W ) and used in equation (5). Consistent 
Lr SS 

with the above manipulations, it is also convenient to define a hold- 
up ratio 



''wj ■ W Lj^ W Lr^SS 

Thus, the model equations (1) and (2) now become 



( 6 ) 



“ *ij = R^T [ (SinpUtS - £° ut P uts > m a SS - Vij ] < 7) 

IT 1 = = ^ [ (Einput - Ioutput) energy ' ^Wj h j] < 8) 



Variation in liquid holdup is allowed for by using expressions 

describing stage hydraulic behavior. Several detailed and very good 

correlations describing such behavior have been presented recently, 

i.e. see papers by Waggoner [16], Tetlow [17] and Peiser and Grover. 

[1] A modified version of the Peiser and Grover equation was adapted 

for use here as it provided an adequate but simple expression of the 

form W * 'p(L.). See Appendix I for derivation. The equations used 
L J J 

are 






nrp. 



Lj 



+ 3 p. 



i J 



WHERE 



<Vss 






( 9 ) 



3 = 



^ W Lr^SS 



cyil ij 
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*Wj = “Plj + 



v ? -Vs,. 



(9a) 



An overall mass balance yields the following additional equation 

EL. . = (Einput - Eoutput) .. .... 

Wj mass overall (10) 



Several additional relations are needed to account for the coef- 
ficients of equations (7) and (8) varying with time if any numeri- 
cal solution of the model equations is to be attempted. These are 
stated below. 

The equilibrium relations are assumed to be of the form 



K ij = exp 



(\ 

T. h B i + c i T j 

J J 



(ID 



Normally is a function of temperature, pressure and composition. 

In equation (11) it is assumed that there is no composition de- 
pendence and pressure drop through the column is negligible. 

The individual component enthalpies are given by 



h i ■ a i T j + b i 



( 12 ) 



where the constants a^ and b^ are specific for each component. This 
same form of temperature dependence is also assumed to hold for 
individual component vapor enthalpies, H^, employing different 
constants. The values for these constants were obtained by cross- 
plotting data from Maxwell. [20] Consistent with the preceding 
relationships for enthalpies, the enthalpy of the liquid on stage 
j, designated h^, is considered a function of composition and 
temperature, viz. 



h. 

J 



f (X, 



V 



(13) 



For other than steady state conditions, X 



ij 



and T. are functions of 
J 
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time so that the derivative of h 



j 



can be written 




( 14 ) 



Assuming that h_. can be expressed as a linear combination of 
individual component enthalpies, namely 



* h . 

3 0 






(15) 



it is possible to evaluate the partial derivatives in equation (14) 
and substitute back giving 



h. 

J 



f dX 1 


i r 


IPiW VJ 


1 + |Pi x ij 




Using equation (3) in the form 



(16) 



PiVj = 1 



(17) 



and taking the time derivative gives 




(18) 



(19) 



Substituting (19) and (18) into (16) and converting to dimensionless 
time, t, gives 
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( 20 ) 



v piVV^j + Pi x ij 



( giAi 



Fi/ij'w • c i> 



The changes in composition from the integration of the mass 
balance equations (7) across At give values of and at 
(t + At). From the X^, the T at (t + At) are obtained by a bubble 
point calculation. Thus, equation (20) provides an independent 

estimate of h. since all of the terms are known at (t + At). 

J 

Equation (20) used in conjunction with equations (8), (9) and (10) 
suggests an iterative scheme to correct internal column flows * 
across time. 

In summary, then, equations (4), (7) through (12), (15), and 
(20) are the model equations used for a numerical solution of the 
general stage dynamic response. 

It might be noted here that the individual component mass balance 
equation* (7) can be considered to be a finite difference approxi- 
mation to (M*N + 2) coupled, nonlinear, partial differential equations 
in X_, the independent variables being time and distance along the 
column. These partial differential equations are of second order in 
distance along the column and first order in time and also contain 
some first order terms in distance. This would indicate that the 
solutions do indeed describe traveling waves of composition changes. 
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3. Numerical Solution. 



Preliminary Discussion 

The data and results of this study were obtained using a CDC 
1604 digital computer to simulate the column. This section will treat 
the numerical integration of the general stage equations. The actual 
programs used for various boundary conditions are given in Appendix 
II. 

Before it is possible to simulate a transient response, it is 
necessary to know the initial steady state conditions for the nu- 
merical integration. All initial and final steady state conditions 
were calculated using Program I given by Hanson et al. [11] Results 
for the unperturbed and perturbed steady state conditions are tabu- 
lated in Appendix VII. Also, since the constants in the hydraulic 
equations (9) contain specific column parameters, some column design 
parameters for a typical column must be specified. See Appendix III 
for values used in this work. 

General Numerical Scheme 

A general approach for obtaining a numerical solution to a system 
of nonlinear differential equations with time varying coefficients is 
to set up a trial and error sequence for the integration with correct- 
ive iteration of the coefficients eventually forcing the equations to 
be satisfied. As was mentioned in the section on the specific model, 
equations (8) and (20) suggest an iterative scheme for correcting 
internal column flows within each time increment. Following this 
approach, the logic sequence used in programming a numerical solution 
for the response is as follows: 

(a) one stage is designated as the starting point for beginning 
the integration. 

(b) equations (7) for that stage are integrated across one At 
using the Runge-Kutta-Gill algorithm as presented by Gill. [10] The 
coefficients are held constant at their last known values. This 
moves the X„ to the point (t + At). 

(c) a bubble point calculation using equations (3), (4) and (11) 
is performed to establish a new T ^ . 

(d) liquid and vapor stream enthalpies are calculated from 
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equations (12) and (15). 

(e) equations (8) and (20) are evaluated. 



(f) the test 



h j eqn(20) ' h j eqn(8)| 



< e (an arbitrary small 



number) is now used for convergence within a single At. If the test 

fails, for the first iteration only, R is incremented by an arbi- 
• W J i 

trary amount, and is evaluated as Z&^/At. is calcu- 

lated from equation (9), V. is calculated from equation (10), and 
steps (b) through (f) are repeated. 

(g) for the second and subsequent iterations, the last two known 
values of are used to predict by the method of "reguli falsi” 

(21) that value of R . which will cause the condition in (f) to be 
2 ” J 

satisfied. When this condition is met, it is possible to proceed 
to the next stage and so on up and down the column until all stages 
have been integrated. 

(h) increment At and begin sequence again at (a). A flow diagram 
of this procedure is given in Appendix IV. 



Boundary Conditions 

In general, some control must be exercised or the column will 
not operate. For initial investigations it was assumed that the feed 
rate is flow controlled and the feed is saturated liquid at its 
bubble point; also, top product, bottom product and reflux stream 
flows were set and held. Thus, any changes in internal column flows 
are absorbed as holdup changes entirely. In a later section when 
control is discussed, these restrictions are modified. But, for 
the present, these assumptions are satisfactory for checking model 
accuracy and looking at the open loop response. 



This was found to be a good approximation after comparing results 
using equation (9a) and using an arbitrary . Use of equation (9a) 
requires a more complex iteration scheme whicn J is very sensitive to 
changes in internal flows and fails to converge except at very small 
increments of At. 



In steps (f) and (g) several other similar schemes were tried 
which involved equating equation (8) to equation (20) and solving 
for L . R^ and . were then calculated directly using the new 
value^of L.. All sdch attempts failed to produce convergence within 
a single ^ At. 
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A total condenser was used in all simulation runs, both with 
and without an accumulator. Thus, reflux liquid is at its saturation 
temperature at all times. When an accumulator was incorporated, it 
was assumed to be an adiabatic mixer undergoing composition, temper- 
ature and holdup changes which in turn specify reflux stream con- 
ditions . 

The reboiler was treated as a large capacity equilibrium stage 
with a heat exchanger. No heat exchanger dynamics were included, 
but can be incorporated at any time. The reboiler duty was treated 
in two ways. 

(a) Reboiler duty (Q ) was set and held at the new known steady 

R 

state value during the entire transient response. 

(b) Qr> was considered a variable and allowed to float within a At 

R 

increment toward a new value which would correspond to the new equi- 
librium conditions in the reboiler at that point in time. 

The first situation was adopted for runs presented here. Physi- 
cally, this means that the controller set point for energy input to 
the reboiler is not changed. 



Delays 

A major use of the model response will be to derive transfer 
functions for control studies using the analog computer. In the 
many possible control schemes, it may well be that errors are detected 
at small values of elapsed time. Any time delay in the fluid trans- 
port at the column ends will affect the derived transfer function, so 
the model should have transport lag effects incorporated in it. One 
point where these lags occur are in the vapor flow from the top stage 
through the condenser-accumulator and reflux return line. Another is 
in the vapor flow in the reboiler return line. 

To introduce delays at these two points, a simple queing tech- 
nique was used. A plug flow situation was assumed in which the 
entering vapor stream variables (Vy T^, Hy y_) are those leaving 
the top plate or reboiler at a particular increment in time and that 
these variables can be moved across the delay interval without mixing. 
Thus, the variables leaving the delay interval correspond to con- 
ditions "n" At's ago. Diagrammatically, this can be represented as 
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follows : 



(Vj T J H j y ij ) 









At o 


A Ti 


At 3 


7 t 


At 

n 


— (v j H j Vl 
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4. Numerical Results 



Column Response to a Feed Perturbation 
To check the model, computer runs were made with a perturbed 
feed and variable end conditions. See Fig. 1 for physical system. 
The feed entering plate 4 consists of: 



Hydrocarbon 



C 

C 

C 

C 



3 

4 

5 

6 



Original 

Composition 

.25 

.25 

.25 

.25 



Perturbed 

Composition 

.25 

.30 

.20 

.25 



Variable boundary condition employed: 

(a) model with total condenser, constant R^, no accumulator 
and no delays.^ 

(b) model as in (a) with variable added. 

(c) model as in (b) with accumulator added. 

(d) model as in (c) with condenser and reboiler delays added. 

2 

Moczek, Otto and Williams [4] reported a substantial effect of 
the computation time increment used on the accuracy of the transient 

response; this was reflected to a large extent in their derived 

3 

approximate transfer functions. Rosenbrock also discussed this 
point earlier. The main idea is the importance of using a compu- 



For a more complete discussion of this model, refer to a paper 
by Duffin, J. H. , ’’Improving the Accuracy of the Dynamic Response 
of a Multistage-multicomponent Column Model," AIChE Symposium on 
Systems and Process Control , 56th National Meeting, San Francisco, 
California, May, 1965. 

2 

Moczek, J. S., Otto, R. E. and Williams, T. J., "Approximation 
Models for the Dynamic Response of Large Distillation Columns," 
AIChE Symposium on Process Control and Applied Mathematics , Vol. 61, 
1965, pp. 136-146. 

3 

Rosenbrock, H. H. , Brit. Chem. Eng. , Vol. 3, 1958, 
pp. 364-367, 432-435, 491-494. 
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pp. 364-367, 432-435, 491-494. 



24 



FIG. 1 



8 STAGE COLUMN, TOTAL CONDENSER, 
EQUILIBRIUM REBOILER 




EQUILIBRIUM REBOILER 





I 

i 



25 




tation interval of a certain minimum size even though computer time 

is correspondingly increased. The interval may be varied over the 
course of a computation to give desired accuracy and speed up solution 
time at larger values of elapsed time. 

For this study, the computation interval (At) was varied from 
.01 to .3 (comparable to the .02 to 3.0 minute interval used by 
Moczek, Otto and Williams) and very little effect on accuracy was 
noticed. This was attributed to the small number of plates used 
in this model. See Fig. 2. In Fig. 2 the short, solid, vertical 
lines indicate the spread of points for At varying between .03 and 
.30; the dotted lines join the data points for At = .1. The value 

of At = 0.3 was found to be the upper limit for a stable solution. 

A At of .1 was selected for use in all subsequent runs although 
computer time to attain steady state was on the order of 45 minutes with 

ample print out data for all runs. No attempt was made to vary At 

during a transient solution. 

Figures 3-8 are temperature, mass transfer, liquid flow and 
holdup response curves representing several stages in the column for 
varying column end conditions. These runs provided a check on the 
accuracy of the model in the sense that the model did move toward a known 
steady state with the mole fraction sums always very close to unity 
and all time derivatives tended to zero within the error limit 
allowed by the iterative scheme. Table 1 lists representative values 
of these accuracy criteria over a typical solution interval. 

The effect of including tray hydraulics, lags and delays should 
be a slowing of the mass transfer response with the final steady state 
being the same for all cases. Referring to Fig. 5, the temperature 
response of the top and bottom stages, and more specifically Fig. 3 
and 4 > top stage and feed stage composition response for the 
fraction, it is evident that this is indeed the case. The responses 
flatten out and the approach to equilibrium becomes slower as the 
model is altered to include additional effects. An enlargement of 

^Computer run time was substantially reduced by reprogramming 
in FORTRAN 63 vice 60. 
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the initial portion of Fig. 3 which is shown in Fig. 9 revealed some 
dead time except for the constant holdup model. With only eight 
stages, dead time was very small and noticeable only three to four 
stages away from the perturbation. For all intents and purposes, 
responses were immediate within the column. 

Referring to Fig. 6 and 7, the liquid flow responses on the top 
stage, feed stage and bottom stage were noticeably different from 
the mass transfer responses. The addition of hydraulics, an accumu- 
lator and then delays accelerated initial changes in the liquid flows. 
At larger values of elapsed t, the response curves crossed so that 
final approach to equilibrium was again slower as more effects were 
added to the model. It can be rationalized that, while the mass 
transfer responses show expected behavior, the lags introduced in the 
condenser-accumulator-reflux cycle and the speed of propagation of 
disturbance waves in liquid-vapor flows resulting from the inclusion 
of hydraulic interactions, cause changes in equilibrium conditions 
which badly upset the vapor-liquid flows initially. The rapid shift 
in internal flows should smooth out as the disturbance waves damp out. 
Also, the iterative scheme for movement of column variables implicit- 
ly alters vapor-liquid flows to meet the convergence test within each 

At. This produces rapid initial changes in flows which gradually 

» * 

smooth out as the time derivatives (h^, R^) pass through their 
maxima and become small. The final approach to equilibrium should be 
similar to the mass transfer responses. 

Representative effects of hydraulic equations on feed stage 
holdup are illustrated in Fig. 8. Table 2 gives the magnitude of 
holdup changes for all eight stages. It should be noted that care 
must be exercised in interpreting responses at very low elapsed t 
values. Any mismatch in the steady and unsteady state model will 
cause initial disturbances which must be allowed to die out before 
any desired disturbances is imposed on the system. 

Column Control Studies 

Control in the digital sense was done by letting the computer 
correct the responses through a control variable. This variable was 
adjusted via a control algorithm based on a generated error. The 

•T'nt l u I. 6 () jt'uLlVc W.-IN l i ' Il.i Lilia .i.1.1 .OutKi C.ilM > > ■ • 'I 

»Vrr.iUM<l < a » > » « > » • • 1 » ■ « . . ■ 
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FIG 3 

MASS TRANSFER RESPONSE OF C 4 FRACTION TO FEED PERTURBATION 
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TOP STAGE LIQUID RESPONSE TO FEED PERTURBATION 
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control objective was to maintain constant composition in the 
overhead product by varying the reflux to compensate for disturbances. 

Accordingly, the top product response to reflux changes was 
needed to determine the proper algebraic form of the control equation. 
Also, these open loop responses are necessary for synthesizing transfer 
functions for analog computer studies. Fig. 10 shows the actual model 
response to reflux perturbations. Fig. 11 is a plot of reflux rate 
versus overhead composition at final steady state. Results of 

these runs indicate that operation at a reflux rate around 1.0 limits 
the control action severely. In the vicinity of a 1.0 reflux rate, 
the steady state composition goes through a maximum. Thus, reflux 
changes would be ineffective in compensating for decreases in overhead 
composition at that point. As long as the perturbations into the 
system produce increases in overhead composition, then control can 
be maintained in this region. By shifting the initial operating point 
away from 1.0 reflux rate, fairly effective control can be realized. 

It was interesting to note also how fast the reflux disturbances 
were propagated through the column. Fig. 12 shows that there is a 
small delay before any change in liquid flow is detected in the bottom 
stage and a slightly longer delay before any composition change occurs. 
This 4s expected. It takes a finite amount of time for disturbance 
waves to travel the column length. If the model is performing properly, 
with only eight stages the delay should be slight, but detectable. 

For this control study the model contained the following: 

(a) variable holdup 

(b) delays in condenser and reboiler lines 

(c) large capacity accumulator 

(d) liquid level control in accumulator and reboiler 

The computer was programmed to use a proportional plus integral 
control equation of the form r = r 0 + KpS + J> 

V?here the error, e, was generated by comparing the instantaneous 
value of in the accumulator to the initial starting value. An 
arbitrary error band of .001 in mole fraction was imposed, i.e. 
control action was taken if |x - ^gsl > X denotes composition 

of in the accumulator. Fig. 13 shows the results of varying K p 
and Kj. The values = 1, gave an underdaraped response with a 
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TABLE 1 

MOLE FRACTION SUM a DERIVATIVE VALUES ON 

FEED STAGE FOR A TYPICAL COMPUTER SOLUTION 

(Ar=O.I) 



MODEL WITH HYDRAULICS MODEL WITH HYDRAULICS a 

' 1 ACCUMULATOR 


#0F 

T*S 


M 

X 


h ' 

EONS 


h 

EON IS 


R w 




X — 

w 


• 

h 

EON e 


h 

EON IS 


^ w 


1 


.99002 


- 1 09.9 


- 190.3 


-3 

1.2X10 




.99632 


-189.9 .. 


| 

-190.3 


i.exio* 3 


10 


1.00000 


- 92.2 


-91.0 


-4 

9.5X10 




10000 


-92.2 


-91.9 


9.7XI0' 4 


00 


.99999 


-35.6 


-35.9 


-4 

3.0X10 




1.00003 


-37.5 


-37.1 


4.9XI0 4 


100 


1.00000 


-24.3 


-23.0 


-4 

5.0X10 




1.00001 


•24.6 


•24.1 


4J0XI0* 4 


500 


.99999 


-3,9 


-3.6 


-4 

I.OXIO 




.99999 


-4.7 


•4.2 


1.8 XI0‘ 4 


1000 


.99999 


-.7 


-.7 


-5 

8.0X10 




1.00000 


-.5 


-.6 


-8 

7.0X10 


♦500 


.90009 


♦.e 


* .2 


-5 

e.oxto 




.99999 


• .4 \ 


-.4 


-5 

o.oxto 



TABLE L 

INITIAL & FINAL VALUES OF HOLDUP RATIO 



STAGE # 


REBOILER 


1 


2 


3 


4 


5 


6 


7 


8 


R w INITIAL 


10.0000 


.9059 


.9946 


.9991 


1.0000 


.8442 


.8641 


.8826 


.899 


R w FINAL 


10.2100 


.9822 


.9971 


1.0087 


1.0145 


.8673 


.8919 


.9050 


.9IIC 


Ar w 


>•2100 


-.0037 


.0025 


.0096 


.0148 


.0231 


*0272 


.0224 


.0111 



40 



long settling time. Increasing the proportional gain to five, 
resulted in an overdamped response which settled out very quickly. 
However, a = K^. = 10 produced an unstable response with increasing 
amplitude of oscillation. No attempts were made to optimize the 
response in any way. This portion of the work was merely to ascertain 
if control could be done in a reasonable time using digital control. 
Results were encouraging in this respect and the model performed 
adequately in all cases. 



5. Transfer Function Synthesis. 



Technique Used 

The purpose of this section is to outline the technique used 
in synthesizing an approximate transfer function for the column and 
present the results of one specific case. 

The Laplace transform of X(t), 

<oo 

X(s) = e‘ St X(t)dt (21) 

J 0 

can be transformed to the interval (0,1) by means of the substitution 



Then 

X(s) = I r S ‘ 1 X(-log r)dr (22) 

*'o 

This integral can be approximated by a Gaussian quadrature formula 
of order N appropriate to the interval (0,1), 



X(s) 




w i X(-log r i )r i S ‘ L 



s 



1 , 2 , 



N 



(23) 



where (r^) are the roots of the shifted Legendre polynomials, 
P N* (X)l 

P N* (X) = P N (1 ‘ 2X) 



and (w^) are the corresponding Chris toffal weights. 

^Tabulated values of the (r.) and the (w.) are given for N up 
to 15 in Bellman', et al. [15]. 1 See also Beliraan, et al. [12] for 
general approach to numerical evaluation of Laplace transforms. 
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Use of equation (2.3) in its present form restricts the X(t) 

values to a very small interval. The upper limit of the interval does 
not increase significantly as the order of N is increased, but the 
computational effort does. Therefore, to expand the X(t) interval by 
a factor, n a," one makes use of the multiplicative property of the 
transform and equation (23) becomes 



X(s/a) 



E w.X(-a log r . )r 
i=l 1 11 



5-1 



(24) 



= 1 , 2 , 



N 



If the functional form of X(s) is known approximately and it is 
desired to improve the estimate of parameters contained in X(s), one 
can use a least squares curve fitting technique on the values of 
X(s) from equation (24) and the calculated values of the approximate 
or assumed functional form of X(s), i.e. 

S = - X < S >Lumed} 



where S is to be minimized with respect to the parameters in 
X(s) 

assumed 

Curve fitting in the "s" domain in lieu of the time domain is 
preferred since the functional form of typical transfer functions 
for many complex processes involving heat and mass transfer have 
been investigated and reported in the literature. See for instance 
references [13], [18] and [19]. An alternative route for getting 
an approximate transfer function is to plot the normalized model 
response on semilog paper. The equation of the best asymptote to 
this curve yields the first terra in the approximate equation of the 

response, i.e., X(t) = a^ t ^ bl + + Next the 

difference between this asymptote and the actual response is plotted 
on semilog paper; the equation of the best asymptote gives the second 
terra of the equation for X(t). In a like manner, the differences 
are plotted until as many terms as needed are obtained for the degree 
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of accuracy desired. A least squares fit could then be used as be- 
fore 
X(s) 



fore to optimize the parameters in X( s ) a pp rox respect to 
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Transfer Function for Top Product Response to a Step in Feed 



The model contained an accumulator, delays in condenser and re- 
boiler, and liquid level controls in reboiler and accumulator. 

The transfer function G(s) for the situation 



STEP CHANGE 
IN FEED 






G(S) 



V s > 



TOP PRODUCT 



was the case studied. G(s) was assumed to be of the form 



“V 

f v Ke 

6CS) = ( Ti 8+1)(t 2 s+1) 

and the values of K, t^, and were estimated graphically as 
described by Moczek et al. [13] from the open loop time responses. 
The results obtained were as follows: 



K 

1 



T d T I t 2 

3 tau 36.7 tau 4.1 tau 

units units units 



1 



.08 min 9.8 min 



1.09 min 



If one denotes F(s) as the perturbation and X(s) as the response, 
the G(s) = X(s)/F(s). Thus, for a step in composition of the 
feed of 0.05, F(s) = .05/s and the functional form of X(s) to be 
used with equation (24) becomes 



X(s) 



nc ‘ • 3S 

.05 e 

s(36. 7s+l) (4. ls+1) 



A factor of 100 was used to expand the useful X(t) interval so 
that the actual form of X(s) used in the least squares fit was 
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X(s/100) 

100 



s = 1,2 



N 



The least squares program used to optimize the three parameters 
(K, t 2 ) is given in Appendix V. The results for a three para- 
meter fit were K = .0536, = 36.2004 and t 2 = 4.4585 or 



G(s) 



1.072 e~ * 3s 

(36. 2004s+l) (4.4585s+l) 



(25) 



Using equation (25), X(s) - .05/sG(s) was inverted back to the time 
domain and plotted along with the model response for X(t). The 
results are shown in Fig. 14. The fit is excellent with only small 
deviations at large values of elapsed time--well beyond the control 
region. 
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10 X 10 to the ^ inch, 5th lines accented. 
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6. Summary and Recommendations. 



A dynamic model has been developed which performs exceedingly 
well under the assumptions made in the mathematical development. 
Additional effects can easily be incorporated into the present model 
once the mathematics are specified. Preliminary control studies 
using the model for error generation indicate single loop and multi- 
loop studies can be done using a reasonable amount of computer time 
and appropriate transfer functions can be obtained from the model 
responses. All computation time can be substantially reduced by 
converting the present Fortran 60 programming to Fortran 63 or 
Fortran 4. Thus, the model can be a useful tool for design and con- 
trol work. 

Future work might be directed toward inclusion of pressure 
effects in the column, downcomer effects, reboiler dynamics and a 
partial condenser to extend the range of control variables and 
interactions which affect the system. Data on an actual operating 
column would be quite helpful in this respect, so that the model 
could be tailored to the actual plant conditions. In this way control 
studies and model responses could be correlated and checked against 
an actual operational plant in hopes of improving the design and 
control loops. 



47 



7. Acknowledgements. 



The author is grateful to Professor J. H. Duffin for proposing 
the research and for his help in completing the project. 



48 



BIBLIOGRAPHY 



1. Peiser, A. M. and Grover, S. S., ’’Dynamic Simulation of a 
Distillation Tower,” Chemical Engineering Progress , v. 58, 

No. 9, pp. 65-70 (1962). 

2. Rosenbrock, H. H. , "The Control of Distillation Columns," Trans . 
Inst, of Chem. Engrs. (London) 40, pp. 35-53 (1962). 

3. Sargent, R. W. H. , "The Dynamic Behavior of Multistage Systems: 
Further Improvements in the Numerical Solution," ibid. , 41, 

pp. 51-60 (1963). 

4. Moczek, J. S., Otto, R. E., Williams, T. J., "Approximation 
Models for the Dynamic Response of Large Distillation Columns," 
Automatic and Remote Control , Butterworth, London, 1963. 

5. Anisimov, I. V., "A Study of the Dynamic and Static Character- 
istics of the Process of Fractional Distillation," ibid. , 
Butterworth, London, 1963. 

6. Voetter, H. , Plant and Process Dynamic Characteristics , Butter- 
worth, London, 1957. 

7. Rijnsdorp, J. E. and Maarleveld, A., "Use of Electrical Analogs 
in the Study of the Dynamic Behavior and Control of Distillation 
Columns," Proc. of the Joint Symp. on Instrumentation and 
Computation in Process Development and Plant Design , London 
(1959). 

8. Williams, J. J., "The Status of Studies of the Dynamics of Mass 
Transfer Operations--A Review and Commentary," Chem. Eng. Prog. 
Symp. , Series #46, v. 59 (1963). 

9. Mah, R. S. H. , Michaelson, S., and Sargent, R. W. H. , " Chemical 
Engineering Science ," 17, pp. 612-639 (1962). 

10. Gill, S., Proc. Cambridge Phil. Soc. , 47, p. 96 (1951). 

11. Hanson, D., Duffin, J. , and Somerville, G., Computation of 
Multistage Separation Processes , Reinholdt, New York, 1962. 

12. Bellman, R., Kagiwada, H., and Kalaba, K. , IEEE Trans, on Auto - 
matic Control , v. AC-10, p. Ill (Jan., 1965). 

13. Moczek, J., Otto, R. , and Williams, T., "Approximation Models 
for Dynamic Response of Large Distillation Columns," Process 
Control and Applied Mathematics , CEP Symposium, Series #55, 
v. 61, p. 136 (1965). 



49 



14. 



Duff in, J. H. , "Improving the Accuracy of the Dynamic Response 
of a Multistage-Multicomponent Column Model," AIChE Symposium on 
Systems and Process Control , 56th National Meeting, San Francisco, 
May, 1965. 

15. Bellman, R.j Kalaba, K., Prestrud, M., and Kagiwada, H. , 

Invariant Imbedding and Time Dependent Transport Processes , 
American Elsevier Publishing Company, Inc., New York, 1964. 

16. Waggoner, R. C., AIChE Journal , v. 11, p. 112 (1965). 

17. Tetlow, N. J., "Analysis and Control of a Generalized Model of 
a Continuous Distillation Column," Doctoral Dissertation, Texas 
Agricultural & Mechanical College, College Station (1966). 

18. Harriot, P., Process Control , McGraw-Hill Book Company, Inc., 

1964. 

19. Williams, T. J., Systems Engineering for the Process Industries , 
pp. 67-69, McGraw-Hill Book Company, Inc., New York, 1961. 

20. Maxwell, J. B., Data Book on Hydrocarbons , D. Von Nostrand 
Company, Inc., New York, 1951. 

21. Lapidus, L., Digital Computation for Chemical Engineers , McGraw- 
Hill Book Company, Inc., New York, 1962. 

22. Rosenbrock, H. H., Brit. Chern. Eng. , v. 3, pp. 364-367, 432-435, 
491-494 (1958). 



50 



APPENDIX I 



DEVELOPMENT OF HYDRAULIC EQUATIONS 



The desired functional form of the expression relating liquid 
holdup to column flows is W = f(Lj). The following development 
is adapted from an analysis for perforated plates done by Peiser and 
Grover [1]. 




Wj ^ 
? T j s 

A tj = 

P Lj" 
h = 

1. = 
J 

8 = 
V = 
L = 



weir height(ft) 

liquid height on the tray(ft) 

liquid height in downcomer (ft) 

2 

cross sectional area of downcomer(ft ) 

2 

cross sectional area of tray(ft ) 

3 

liquid molar density(raoles/f t ) 
height of the weir crest(ft) 
width of the weir(ft) 

2 

gravitational constant (ft/sec ) 
molar vapor f low(moles/time) 
molar liquid f low(raoles/time) 
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s 

°Dj = 
C Lj = 
\d = 

v . = 



D 

«v s 



correlation constant in the Francis weir equation 

liquid flow resistance in downcomer 

liquid flow reversal resistance in downcomer 

pressure drop for vapor flow and liquid head on tray abobe 

liquid velocity in downcomer 

locally defined constant 

area under downcomer 

vapor velocity at tray (j + 1) 

perforation area 

total number of perforations at tray (j + 1) 
orifice discharging coefficient (approx. 0.6) 
molecular weight of liquid 
molecular weight of vapor 



Using the Francis weir formula for large straight weirs 

Va. 

L j = kh (1) 



where 



cA? i .P Lj 



From the diagram 



h = 9 Tj - fyj 



( 2 ) 



Using a pressure balance 



= + + 'ppj 



( 3 ) 

( 4 ) 
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